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A  TECHNICAL  DESCRIPTION  OF  THE 
NAVDAS  ADJOINT  SYSTEM 


1  Introduction 

Adjoint  versions  of  atmospheric  forecast  models  have  been  developed  and  applied, 
beginning  in  the  mid  1980s,  to  a  wide  variety  of  applications  in  numerical  weather 
prediction,  including  data  assimilation,  sensitivity  studies,  singular  vector  generation 
and  targeted  (or  adaptive)  observing.  The  literature  describing  adjoint  methods  in 
numerical  weather  prediction  is  now  quite  extensive,  and  a  few  selected  papers  on 
these  topics  include  Talagrand  and  Courtier  (1987),  Buizza  et  al.  (1993),  Langland 
et  al.  (1995),  Rabier  et  al  (1996),  Errico  (1997),  Gelaro  et  al.  (1998),  Palmer  et 
al.  (1998),  Bergot  et  al.  (1999),  Langland  et  al.  (2002),  Bergot  and  Doerenbecher 
(2002),  and  Leutbecher  et  al.  (2002).  More  recently,  it  has  been  shown  (Baker  and 
Daley  2000)  that  the  adjoint  of  a  data  assimilation  procedure  can  be  used  to  provide 
the  sensitivity  of  a  forecast  costfunction  to  the  observations  and  background  in  an 
atmospheric  analysis.  This  extension  of  adjoint  sensitivity  methods  into  observation 
space  is  a  significant  advance  in  efforts  to  understand  factors  that  control  atmospheric 
predictability  and  to  improve  data  assimilation  procedures. 

The  first  version  of  the  NAVDAS  (NRL  Atmospheric  Variational  Data  Assimila¬ 
tion  System)  adjoint  was  developed  by  Roger  Daley  during  1999-2000  and  is  briefly 
described  in  Daley  and  Barker  (2001).  This  code  was  used  for  the  Ph.D.  research 
of  Nancy  Baker  (Baker  2000,  Baker  and  Daley  2000)  and  for  preliminary  observa¬ 
tion  sensitivity  tests  by  Rolf  Langland.  However,  substantial  revisions  to  the  regular 
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(forward)  NAVDAS  code  between  2000  and  2002  made  it  necessary  to  develop  a  new 
version  of  the  adjoint,  and  a  code  development  project  was  therefore  initiated  in  2002. 
The  current  NAVDAS  adjoint  was  primarily  developed  from  a  version  of  NAVDAS 
obtained  on  June  18,  2002.  It  is  a  completely  new  adjoint  code,  although  based  on  the 
framework  of  the  original  Daley  version.  Since  NAVDAS  is  a  linear  procedure,  there 
is  no  separate  tangent  linear  version  of  NAVDAS,  and  the  adjoint  was  developed  line 
by  hne  directly  from  the  regular  NAVDAS  Fortran  code.  The  NAVDAS  adjoint  code 
was  further  updated  by  the  authors  using  NAVDAS  code  obtained  in  March  2003  and 
January  2004  that  included  modifications  to  background  error  covariance,  radiance 
assimilation,  and  other  changes.  This  report  provides  the  first  complete  description 
of  the  NAVDAS  adjoint,  including  observation  sensitivity  equations  and  overview  of 
the  Fortran  code  (Section  2),  steps  involved  in  rimning  and  verifying  the  adjoint  code 
(Section  3),  and  examples  of  observation  sensitivity  applications  (Section  4). 

2  NAVDAS  Adjoint  Development 
2.1  Observation  Sensitivity 

The  basic  purpose  of  the  NAVDAS  adjoint  is  to  provide  the  sensitivity  of  a  user- 

defined  scalar  costfunction^  J  with  respect  to  the  vector  of  observations  y  used  by 
can  be  any  differentiable  function  of  the  forecast  model  variables,  such  as  forecast  error, 
vorticity,  surface  pressure,  temperature,  wind  speed,  etc. 
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NAVDAS  in  its  assimilation  procedure.  The  costfunction  J  is  a  function  of  a  model 
forecast  (x/)  started  from  an  analyzed  initial  condition  (xa)  and  J  may  also  be 
referred  to  as  a  ’’forecast  response  function”  or  ’’forecast  aspect.”  The  ’’observation 
sensitivity”  vector  dJ  /  dy  Ss  &  gradient  in  ’’observation  space”,  e.g.,  its  elements 
exist  at  the  locations  of  observations,  and  it  is  defined  by 


dJ  _  X 


(1) 


where  dJ  /  dxa  is  an  analysis  (initial  condition)  sensitivity  gradient  provided  by  the 
adjoint  of  a  forecast  model  such  as  NOGAPS  (Rosmond  1997,  Hogan  et  al.  1999). 
The  adjoint  operator  is  the  transpose  of  the  Kalman 

gain  operator,  defined  by  K  =  PaH"^  [HP^H^  4-  R]“^  that  represents  the  regular 
(forward)  NAVDAS  procedure  (Daley  and  Barker  2001).  Note  that  [HPfcH'^  +  R]“^ 
is  self-adjoint.  H  is  a  linearized  operator  used  to  project  from  gridpoint  to  observation 
space,  and  P^  and  R  are  the  background  and  observation  error  covariance  matrices, 
respectively.  The  main  components  of  the  observation  sensitivity  calculation  are 
outlined  in  Fig.  1 
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Figure  1:  Components  of  the  observation  sensitivity  calculation  using 
adjoint  versions  of  NOGAPS  and  NAVDAS,  where  is  the  tangent 
propagator  of  NOGAPS,  xy  is  the  forecast  state  vector  and  x^  is  the 
analysis  state  vector.  Other  symbols  are  defined  in  Section  2. 1 . 


The  initial  condition  sensitivity  is  defined  by 


dJ  ^  T  dJ 
9Xa  (9xy  ’ 


(2) 


where  is  the  adjoint  (transpose)  of  the  tangent  linear  propagator  L  of  the  forecast 
model.  Since  the  Kalman  gain  operator  in  NAVDAS  is  hnear,  the  calculation  of 
sensitivity  using  (1)  is  exact.  However,  calculation  of  the  analysis  sensitivity  dJ  /  5xq 
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by  the  NOGAPS  adjoint  (2)  involves  tangent  linear  approximations  and  this  limits 
the  useful  quantitative  application  of  observation  sensitivity  to  short-range  forecasts 
of  approximately  three  days  or  less.  In  addition,  the  current  versions  of  the  NOGAPS 
and  NAVDAS  adjoints  do  not  include  consideration  of  moist  processes  or  moisture 
observations,  which  can  be  an  important  limitation  in  some  applications. 

A  perturbation  of  the  initial  conditions  Sxa  can  be  propagated  to  the  forecast 
time  by 

SXf  =  LSxa  .  (3) 

The  introduction  of  Sxa  changes  the  scalar  costfunction  by  an  amount  SJ,  which 
is  a  function  of  Sxf 


J-t-SJ  =  f  (xf  +  Sxf) 

(4) 

SJ  =  f  (Sxf)  =  f  (L^Xa)  . 

(5) 

Since  SJ  is  a  function  of  Sxa,  we  can  extend  the  definition  of  the  sensitivity 
gradient  in  terms  of  the  observation  and  background  used  by  the  data  assimilation 
procedure.  The  basic  linear  form  of  the  NAVDAS  analysis  can  be  written  as 

Xa  =  Xfc  -)-  K  (y  -  Hxb)  ,  (6) 

where  Xb  is  a  background  forecast.  In  the  regular  form  of  NAVDAS  (Daley  and 
Barker,  2001)  the  operator  H  may  be  nonlinear  for  some  t5q)es  of  observations.  If  we 


5 


perturb  the  observations  by  Sy,  then 


Xa  +  (5xa  =  Xi,  +  K  (y+5y  -  Hxj)  .  (7) 

The  adjoint  of  (7),  which  defines  the  sensitivity  to  observations,  is  (1).  However, 
if  we  instead  perturb  the  background  by  ^x^,  then 

Xa  +  Sxa  =  {xb  +  5xb)  +  K  (y  -  Hxfc  -  H5x6)  .  (8) 


The  adjoint  of  (8),  which  defines  the  sensitivity  to  the  background,  has  two  terms 


dJ  _  dJ 

dXb  dXa 


In  (9),  dJ  /  dxa  is  the  observation  sensitivity.  Therefore  we  may  re-write  (9) 


dJ  _  dJ 

dXb  dXa 


which  is  an  equation  for  the  sensitivity  of  J  to  the  background  in  gridpoint  space.  If 
there  are  no  observations,  it  is  clear  that  dJ  /  dxb  =  dJ  /  dx^  because  the  analysis  Xq 
is  then  simply  the  background  xj  carried  forward  with  no  change  and  dJ  /  5y  =  0. 
The  operator  H"^,  which  is  required  to  project  dJ  /  dy  from  observation  space  to 
gridpoint  space,  is  not  available  as  a  stand-alone  component  of  the  regular  NAVDAS 
code  but  it  is  currently  being  developed  as  part  of  the  adjoint  system. 

Observation  sensitivity  can  be  used  in  a  variety  of  ways  to  study  the  impact  of 
observations  on  forecast  outcome.  For  example,  a  Taylor  series  expansion  can  be  used 


6 


to  estimate  the  effect  of  an  observation  pertiirbation  vector  in  terms  ol  5  J 


/.^  ,2  lavx 


(11) 


If  the  perturbations  represented  by  5y  do  not  exceed  the  magnitude  of  t3T)ical 
innovations,  the  Taylor  series  can  be  tnmcated  to  just  the  l®*-order  term 


The  scalar  quantity  5J  can  represent  the  forecast  impact  caused  by  chstnges  to 
real  observations,  or  even  the  addition  of  hypothetical  observations.  It  will  be  shown 
in  Section  4.2  that  observation  sensitivity  can  also  be  used  to  quantify  the  impact  of 
innovations  (observation  -  background)  on  short-range  forecast  error  differences.  In 
fact,  it  is  valid  to  interpret  dJ  /  dy  as  &  sensitivity  to  the  innovation  vector,  and  Sy 
as  a  perturbation  or  component  of  the  innovation  vector.  Following  (10)  the  impact 
of  backgroimd  perturbations  can  be  evaluated  using  the  expression 

where  the  first  term  on  the  RHS  of  (13)  is  evaluated  in  gridpoint-space  and  the 
second  RHS  term  is  in  observation  space. 
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Figure  2:  Primary  subroutines  in  forward  version  of 
NAVDAS.  Symbols  are  defined  in  Section  2.1. 


2.2  Adjoint  Code  Description 

The  NAVDAS  adjoint  code  is  written  in  Fortran  90  and  designed  to  run  on  multiple 
processors  using  MPI  (message-passing-interface).  The  primary  difference  between 
the  forward  and  adjoint  codes  for  NAVDAS  is  the  transposition  of  the  post-multiplier. 
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Whereas  in  the  forward  NAVDAS  (Fig.  2)  the  post-multiplier  is  appHed  after 
the  solver  [HP^H’’  -I-  R]“\  in  the  adjoint  of  NAVDAS  (Fig.  3)  the  transpose  of  the 
post-multiplier  HP^  is  invoked  first,  and  then  the  solver  (which  is  the  same  in  the 
forward  and  adjoint  codes)  is  applied.  In  addition,  the  NAVDAS  adjoint  uses  special 
routines  for  interpolation  of  the  input  analysis  sensitivity  vector  dJ  /  5xa  and  for 
post-processing  of  the  observation  sensitivity  output. 

The  computational  time  required  to  run  the  NAVDAS  adjoint  is  less  than  that  of 
the  regular  NAVDAS  analysis,  because  the  adjoint  does  not  require  that  the  innova¬ 
tion  processing  of  the  regular  NAVDAS  be  repeated.  The  adjoint  code  is  structured  so 
that  major  functions  are  controlled  by  a  driver  program  and  several  primary  subrou¬ 
tines,  which  are  described  below.  However,  all  details  of  the  adjoint  code  development 
are  not  included  here,  and  the  actual  Fortran  code  should  be  examined  (see  Appendix 
A)  if  additional  information  is  required. 
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Figure  3:  Primary  subroutines  in  adjoint  vasicm  of  NAVDAS.  Symbols 
are  defined  in  Secti(Mi  2.1.  Hie  caret  (^)rqjresents  37/5  (  ). 

2.2.1  adj_nav3d_ nogaps  The  main  program  module  of  the  adjoint  code  is  sim¬ 
ilar  to  the  main  program  of  the  regular  (forward)  NAVDAS  code  (nav3d  nogaps), 
except  that  a  call  is  made  to  navdas_adj,  rather  than  to  navdas.  Certain  para¬ 
meters  and  switches  are  set  differently  for  running  the  adjoint  code.  For  example, 
the  adjoint  does  not  repeat  any  of  the  quality  control  checks  that  are  performed  in 
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the  regular  NAVDAS  run.  The  quality  control  flag  that  results  from  the  innovation 
and  buddy  check  quality  control  is  read  from  the  rsd_vec  file  that  has  been  pro¬ 
duced  in  the  regular  NAVDAS  nm.  The  adjoint  also  reads  a  quality  control  flag 
from  the  innovation  file  that  results  from  quality  control  checks  in  the  NAVDAS 
pre-processor. 

Switch  Settings  in  NAVDAS  Adjoint: 

•  innov_check  =  F 

•  buddy_check  =F 

•  ghost  =  F  (this  is  an  option  for  ensuring  that  preconditioner  prisms  have  suffi¬ 
cient  observations  in  the  regular  NAVDAS) 

•  isentropic  =  F  (the  adjoint  does  not  currently  allow  this  option) 

•  jmin  =  F  (the  jmin  diagnostics  are  not  computed  in  the  adjoint) 

•  analysis_error_exe  =  F  (the  call  to  navdas_aeiT  for  analysis  error  estimates 
is  not  part  of  the  adjoint  code) 

2.2.2  navdas_adj  This  subroutine  reads  the  various  input  files  required  to  run 
the  NAVDAS  adjoint,  calls  the  major  subroutines  involved  in  the  sensitivity  calcula¬ 
tions  and  writes  the  observation  sensitivity  to  an  output  file.  The  first  step  is  to  read 


11 


the  innovation  file  and  rsd_vec  file  produced  by  the  regular  NAVDAS  assimila¬ 
tion.  These  files  contain  information  that  specifies  observation  value,  location,  time, 
background  value,  errors  assigned  to  the  observation  and  background,  quality  control 
flags  and  other  information  valid  at  the  assimilation  time  for  which  the  sensitivity  is 
to  be  calculated.  The  observations  considered  by  the  adjoint  can  be  real  observations 
provided  in  the  innovation  file,  or  they  can  be  hypothetical  observations  defined  by 
the  user. 


Subroutine 

Array(s) 

Remarks 

navdas_adj 

tempsen,  usen, 
vsen 

Analysis  sensitivity  (temperature, 
u-wind,  v-wind)  on  the  1  60-leveI 
“output”  grid  of  NAVDAS. 

adjoint_correction 

adjoint_eigen_grid 

Analysis  sensitivity  in  vertically- 
decomposed  eigenvector  space  (60 
levels). 

adjoint_n3dvar 

eob_work 

Observation  sensitivity  in 
eigenvector  space  after  call  to 
adjoint_post_multiply 

cob 

Observation  sensitivity  in  physical 
space  before  call  to  solver 

navdas_adj 

ob_sens 

Observation  sensitivity  vector  in 
observation  space  after  call  to 
adjoint_n3dvar. 

Table  1:  Arrays  that  contain  analysis  and  observation 
sensitivity  in  the  NAVDAS  adjoint  code. 


The  analysis  sensitivity  vector  dJ  /  dxa  provided  by  the  NOGAPS  adjoint  is  read 
from  a  saved  file  and  put  into  the  arrays  tempsen,  usen  and  vsen,  which  contain 
the  sensitivities  for  temperature,  u-wind,  and  v-wind,  respectively.  The  NOGAPS 
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adjoint  also  produces  a  sensitivity  to  terrain  pressure,  but  this  component  is  not  used 
as  input  for  the  NAVDAS  adjoint,  because  the  forward  NAVDAS  does  not  produce  a 
pressure  analysis  increment  as  an  output  field.  An  array  phisen  appears  in  the  code 
but  would  be  used  only  if  the  NOGAPS  adjoint  produced  sensitivity  to  height  instead 
of  temperature.  Note  that  a  special  interpolation  procedme  (see  Section  3.1)  is  used 
to  interpolate  dJ  /  dxa  from  the  NOGAPS  grid  to  the  NAVDAS  grid  before  running 
the  NAVDAS  adjoint. 

A  call  to  adjoint_correction  converts  the  sensitivity  from  the  analysis  gridpoint 
space  to  eigenvector  space,  and  the  sensitivity  is  then  contained  in  the  array  ad- 
joint_eigen_grid  (Table  1).  The  subroutine  windrotate_grid  adjusts  the  analysis 
sensitivity  winds  from  grid  orientation  to  spherical  coordinates,  and  adjoint  _n3dvar 
is  then  called  to  perform  the  major  part  of  the  observation  sensitivity  calculations. 
Finally,  windrotate_ob  is  called  to  convert  the  sensitivity  winds  (now  in  obser¬ 
vation  space)  back  to  grid  orientation  firom  spherical  coordinates.  The  observation 
sensitivity  is  then  written  to  an  output  file. 

2.2.3  adjoint_correction  This  subroutine  is  called  from  navdas_adj  and  con¬ 
verts  the  analysis  sensitivity  vector  (contained  in  the  arrays  tempsen,  men,  vsen) 
from  NAVDAS  gridpoint  space  to  eigenvector  space  {adjoint_eigen_grid)  near  the 
beginning  of  the  observation  sensitivity  calculations.  The  ’’zero-order”  vertical  in¬ 
terpolation  option  of  NAVDAS  is  used  here,  because  the  60  pressure  levels  used  for 
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the  arrays  tempsen,  usen,  vsen  have  been  specified  to  exactly  match  the  60  levels  of 
adjoint_eigen_gri±  The  pressure  levels  for  temperature  and  wind  are  defined  by  the 
arrays  prestout  and  pressout,  respectively.  Also  in  this  subroutine  the  temperature 
and  wind  sensitivities  are  ’’normalized”  by  multiplying  with  the  square  root  of  the 
background  error  variance  (the  array  rmsvar). 

2.2.4  adjoint^nSdvar  This  subroutine  controls  the  two  major  components  of 
the  observation  sensitivity  calculations,  which  are  the  adjoint  (transpose)  of  the  NAV- 
DAS  post-multiplication  step  and  the  iterative  application  of  the  solver.  At  the  begin¬ 
ning  of  this  subroutine,  the  sensitivity  is  contained  in  the  array  adjoint^  eigen_  grid  on 
the  NAVDAS  one-degree  grid  and  in  vertical  eigenvector  space  using  60  vertical  levels. 
The  first  major  step  is  to  apply  the  operator  HP^  (the  transpose  of  the  forward  NAV¬ 
DAS  post-multiplier  P^H'^)  to  the  analysis  sensitivity  as  a  ’’pre-multiplier”  by  calling 
the  subroutine  adjoint  _post_multiply.  The  call  to  adjoint  post  multiply  is 
actually  made  within  a  set  of  do-loops  indexed  over  the  number  of  analysis  grid 
volumes  and  observation  prisms. 

After  calling  adjoint  _post_^multiply  the  sensitivity  is  contained  in  the  array 

eob_work  (Table  1).  The  subroutine  left _ operator  converts  the  sensitivity  from 

the  eigenvector  space  of  eob^work  into  the  array  co6,  which  is  the  sensitivity  in 
physical  (observation)  space.  The  first  and  second  preconditioners  are  applied  and 
subroutine  genince  is  then  called  to  perform  the  iterative  conjugate  gradient  solver 
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I 

component  of  the  NAVDAS  adjoint,  represented  by  +  R]~^.  The  solver 

subroutine  (genince)  is  the  same^  in  the  forward  and  adjoint  versions  of  NAVDAS, 
but  the  argument  list  in  the  call  to  genince  is  altered,  as  follows: 

•  call  genince  (w,  xiv_ob,  cob, . )  in  forward  NAVDAS 

•  call  genince  (w,  cob,  ob_sens,  ....)  in  adjoint  NAVDAS 

In  the  forward  NAVDAS,  the  solver  input  is  the  array  xiv_  ob  and  it  returns  the 
array  co6,  which  is  in  observation  space  prior  to  the  appUcation  of  right_operator. 
In  the  adjoint  NAVDAS,  the  solver  input  is  the  array  cob  and  it  returns  the  array 
ob_  sens,  which  is  the  observation  sensitivity  in  observation  space.  The  final  step  in 
adjoint _n3dvar  is  a  ’’de-normalization”  of  the  sensitivity,  accomplished  by  dividing 
ob_sens  by  the  array  rmsvar.  The  array  ob_sens  is  then  returned  to  navdas_adj 
for  output. 

2.2.5  adjoint_post_multiply  This  subroutine  is  called  from  adjoint _n3dvar 
and  is  used  to  apply  the  HP^  operator  to  the  sensitivity  vector  in  the  first  step  of 
the  NAVDAS  adjoint  calculations.  It  works  in  eigenvector  space  and  represents  the 
transpose  of  the  operator  PftH^  performed  by  subroutine  post_muItiply  in  the 
forward  NAVDAS.  In  the  forward  NAVDAS,  the  input  for  post_multiply  is  the 

^The  convergence  criteria  used  in  the  descent  algorithm  of  the  solver  are  identical  in  the  regular 
NAVDAS  and  adjoint  NAVDAS. 
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array  eob_work  (the  array  cob  projected  into  vertical  eigenvector  space)  and  the 
output  is  an  array  containing  what  is  essentially  the  analysis  correction  vector  in 
eigenvector  space.  In  the  adjoint  NAVDAS,  the  input  for  adjoint_post_multiply 
is  the  sensitivity  array  adjoint_eigen_grid  and  the  output  is  the  adjoint  version  of 
the  array  €ob_  work. 

2.2.6  Arrays  vised  in  NAVDAS  adjoint 

•  tempsen,  usen,  vsen  -  sensitivity  to  temperature  and  winds  (components  of 
dJ  /  axa)  on  the  1°  (lat-lon)  60-level  NAVDAS  output  grid 

•  adjoint_eigen_grid  -  sensitivity  to  temperature  and  winds  {dJ  /  &Ka)  on  the 
1°  (lat-lon)  grid  decomposed  in  vertical  eigenvector  space  (60  levels) 

•  ob_sens  -  sensitivity  to  observations  {dJ  /  dy)  in  observation  space,  cmrently 
includes  temperature,  wind,  height,  and  brightness  temperature  observations 

•  oberr_sens  -  sensitivity  to  specified  observation  error  in  observation  space,  see 
Section  4.5  for  definition 

•  bkerr_$ens  -  sensitivity  to  specified  background  error  in  gridpoint  space,  cur¬ 
rently  not  implemented 

•  prestout,  pressout  -  defines  the  60  pressure  levels  used  for  input  of  temperature 
and  wind  analysis  sensitivity,  respectively 
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•  xiv_  ob  -  innovation  values  read  from  the  innovation  file  produced  by  NAVDAS 

•  rsd_  val  -  residual  (analysis  -  observation)  values  read  from  the  rsd_vec  file 
produced  by  NAVDAS 

•  num_reject  -  observation  reject  flag  values  read  from  the  rsd_vec  file  pro¬ 
duced  by  NAVDAS 

•  cob_fwd  -  ”cob”  vector  from  the  regular  (forward)  NAVDAS  read  from  the 
rsd_vec  file 

3  Running  the  NAVDAS  Adjoint 

The  NAVDAS  adjoint  code  exists  as  a  single  file  written  in  Fortran  90  called  NAV- 
DAS_adj.f  that  contains  the  main  program  and  all  required  subroutines  and  defined 
functions.  The  source  code  is  compiled  using  a  makefile  that  produces  the  executable 
NAVDAS_adj.exe.  There  is  no  attached  code  library.  The  adjoint  calculation  uses 
an  outer  c-shell  script  called  run_navdas_adjoint.  The  c-shell  script  controls  in¬ 
put  and  output  of  certain  files  used  by  the  adjoint  and  invokes  a  k-shell  script  called 
NAVDAS_adj.ksh  that  invokes  the  NAVDAS  adjoint  executable  file.  The  NAV¬ 
DAS  adjoint  source  code  and  scripts  are  available  in  a  directory  on  the  NRL  computer 
hadley  (see  Appendix  A).  The  procedure  for  calculating  observation  sensitivity  with 
the  NAVDAS  adjoint  includes  the  following  steps: 
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Procedure  for  Calculating  Observation  Sensitivity: 

1.  Run  the  regular  NAVDAS  assimilation  to  produce  analyzed  initial  conditions. 
Save  the  innovation  file  and  rsd_vec  file 

2.  Run  the  full-physics  nonlinear  forecast  model  (NOGAPS)  and  save  the  forecast 
trajectory 

3.  Define  a  forecast  costfunction  (J) 

4.  Run  the  NOGAPS  adjoint  to  produce  the  analysis  (initial  condition)  sensitivity 
dJ  /  &x.a 

5.  Interpolate  dJ  /  5xa  from  the  NOGAPS  gaussian  grid  to  the  NAVDAS  output 
grid  (see  Section  3.1,  below) 

6.  Run  the  NAVDAS  adjoint  to  produce  the  observation  sensitivity  dJ  /  dy 

7.  Run  graphics  or  post-processing  to  plot  observation  sensitivity,  observation 
impact  or  other  products 

3.1  Pre-Processing  (Interpolation  of  dJ  /  dxa) 

The  analysis  sensitivity  vector  dJ  /  9xa  comprises  the  ’’input”  to  the  NAVDAS 
adjoint  and  is  provided  on  a  1°  lat-lon  grid  with  60  vertical  levels  that  correspond 
exactly  to  the  levels  used  in  the  vertical  eigenvector  decomposition  at  the  end  of  the 
regular  NAVDAS  analysis.  This  allows  the  zero-order  interpolation  option  to  be  used 
when  the  arrays  that  contain  dJ  /  {tempsen,  v,sen,  vsen)  are  transferred  into 
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adjoint_eigen_gnd  by  the  subroutine  adjoint_correction. 


W, 


Ap  =  Imb 


Figure  4:  Interpolation  of  an  adjoint  sensitivity  gradient  between 
grids  with  different  vertical  configurations.  Used  for  transfer  of 
analysis  sensitivity  from  NOGAPS  sigma  levels  to  NAVDAS  ou^ut 
pressure  levels  in  the  pre-processing  step  for  NAVDAS  adjoint.  W  is 
an  array  that  contains  the  analysis  sensitivity  values  normalized  to 
represent  pressure  intervals  of  1  mb. 


Before  running  the  NAVDAS  adjoint,  a  special  pre-processing  interpolation  pro¬ 
cedure  is  required  to  transfer  dJ  /  5xa  from  the  NOGAPS  gaussian  grid  and  sigma 
surfaces  onto  the  NAVDAS  grid.  Because  the  adjoint  sensitivity  fields  are  gradients, 
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rather  than  conventional  temperature  or  wind  variables,  this  interpolation  step  must 
carefully  account  for  differences  in  horizontal  and  vertical  resolution  of  the  NOGAPS 
and  NAVDAS  grids.  It  is  critical  that  this  step  be  performed  correctly,  so  that  the 
global  sum  of  the  analysis  sensitivity  is  identical  on  both  grids.  Other  procedures, 
such  as  bi-linear  interpolation,  should  not  be  used  to  interpolate  adjoint  sensitivity 
between  grids  with  different  horizontal  or  vertical  resolutions,  as  they  are  likely  to  in¬ 
troduce  serious  distortion  of  the  sensitivity  gradients.  A  version  of  the  pre-processing 
sensitivity  interpolation  procedure  for  NAVDAS  is  provided  as  a  Fortran  90  file  called 
gradtp.f  (see  Appendix  A). 

The  first  step  in  the  interpolation  of  the  analysis  sensitivity  (temperature  and 
winds)  is  a  transfer  of  dJ  /  dxa  from  NOGAPS  sigma-levels  to  the  60  constant- 
pressure  levels  used  by  NAVDAS,  at  every  point  on  the  NOGAPS  gaussian  grid.  The 
procedure  is  illustrated  in  Fig.  4. 

1.  The  sensitivity  {dJ  /  dxa)  values  on  each  NOGAPS  sigma  level  are  normalized 
by  the  thickness  of  the  pressure  layer  (in  mb)  that  each  level  represents.  Thus,  for 
some  level,  n  with  a  pressure  thickness  (Ap)^,  we  let  Sn  =  {dJ  /  dxa)^  /  (Ap)^. 

2.  The  values  of  of  5^,  for  all  n  levels  of  NOGAPS,  are  used  to  define  the 
elements  of  an  array,  W ,  that  contains  the  normalized  sensitivity  values  at  1  mb 
vertical  intervals  from  1  mb  to  1200  mb.  Thus,  if  sigma  level  n  represents  pressure 
firom  500  -  600  mb  at  a  NOGAPS  gridpoint,  the  value  Sn  is  assigned  to  each  of  the 
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k  elements  of  W  from  A:=500  to  600  mb.  Any  elements  of  W  below  the  NOGAPS 
surfaxje  terrain  pressure  are  assigned  a  value  of  zero. 


Latitude 


Longitude 


■  NOGAPS  T79  grid  O  NAVDAS  1»  grid 


Figure  5:  Example  of  interpolations  performed  on  constant  pressure  surfaces  between  the 
NAVDAS  1®  grid  and  the  NOGAPS  T79  grid  (shown  with  a  resolution  of  exactly  1.5®  here  to 
simplify  the  diagram).  The  weighting  factors  on  arrows  directed  towards  NOGAPS  points  represent 
interpolation  of  conventional  fields  and  add  to  1.0.  The  weighting  factors  on  arrows  directed 
towards  NAVDAS  points  represent  interpolation  of  sensitivity  gradient  fields  and  add  to  0.4444, 
which  is  the  ratio  of  the  areas  represented  by  grid  points  on  the  1°  grid  and  the  1.5°  grid,  0.4444  = 
((1.0*1.0)  ^  (1.5*1.5)).  NAVDAS  1®  gridpoints  represent  smaller  areas  than  NOGAPS  T79 
gridpoints,  and  the  sensitivity  magnitude  is  reduced  through  this  interpolation  to  account  for  the 
difference.  Interpolation  weighting  factors  are  determined  by  procedures  illustrated  in  Figs.  6  and  7. 
In  practice,  all  grid  points  are  used  in  the  interpolations. 
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3.  The  normalized  sensitivity  is  then  transferred  from  the  array  W  onto  the 
60-level  NAVDAS  grid,  in  the  reverse  sense  of  Steps  1  and  2.  Each  of  the  60  levels 
on  the  NAVDAS  grid  represents  a  layer  pressure  thickness.  Thus,  if  pressure  level  m 
represents  pressure  from  520  -  560  mb  on  the  NAVDAS  grid,  the  value  {dJ  /  dxa)^ 
is  the  sum  of  the  k  elements  of  W  from  k=520  to  560.  Note  that  the  pressure  levels 
for  temperature  are  not  the  same  as  for  wind. 

4.  After  completion  of  this  step,  the  values  of  dJ  /  dx^  in  any  vertical  column 
can  be  summed  over  all  sigma  levels  on  the  NOGAPS  grid  and  over  all  pressure  levels 
on  the  NAVDAS  grid  and  the  two  sums  should  be  identical. 

The  second  interpolation  step  is  a  transfer  of  dJ  /  dx^  from  the  gaussian  grid 
of  NOGAPS  to  the  1°  lat-lon  NAVDAS  grid,  on  each  of  the  60  levels  defined  in  the 
previous  step.  This  horizontal  interpolation  is  designed  to  conserve  the  global  sum  of 
the  sensitivity  gradient  when  it  is  transferred  between  grids  with  different  resolutions. 
The  ’’forward”  interpolation  is  a  transfer  of  conventional  fields  (temperature,  wind) 
from  the  NAVDAS  grid  to  the  NOGAPS  grid.  As  shown  in  Fig.  5,  the  weighting 
factors  on  arrows  directed  from  NAVDAS  to  NOGAPS  points  (in  the  forward  inter¬ 
polation  sense)  add  to  1.0,  which  conserves  the  values  of  the  conventional  variables. 
In  contrast,  the  weighting  factors  on  arrows  directed  from  NOGAPS  to  NAVDAS 
points  (in  the  adjoint  interpolation  sense)  add  to  0.4444,  which  is  the  ratio  of  the 
areas  represented  by  grid  points  on  the  NAVDAS  1°  grid  vs.  the  NOGAPS  T79  grid 
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('"1.5°).  Thus,  the  sensitivity  magnitude  at  individual  grid  points  on  the  1°  grid  is 
less  than  on  the  T79  grid,  because  each  1°  gridpoint  represents  a  smaller  area,  but 
the  global  sum  of  the  sensitivity  over  all  points  on  each  grid  is  identical.  Different 
weighting  factors  would  be  used  when  the  resolution  of  either  the  NOGAPS  gaussian 
grid  or  the  NAVDAS  grid  is  changed. 


Latitude 


28.5N 


The  distance  from  27. 5N  to  28.45N  is 
assigned  to  the  NAVDAS  grid  point  at 

28. ON.  Its  weighting  factor  is:  (28.45- 
27.5)^  1.50  =  0.6333. 

i 

27.5N 

i 
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The  distance  from  26.95N  to  27.5N  is 

assigned  to  the  NAVDAS  grid  point  at 

27.0N.  Its  weighting  factor  is:  (27.5- 
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r 

26.95)^-1.50  =  0.3667. 

26.5N 

□  29.2N 
O  29.0N 

-  28.45N 

O  28.0N 


4 - O  27.0N 

-  26.95N 

□  26,2N 
O  26.0N 


■  NOGAPS  T79  grid  O  NAVDAS  1®  grid 


Figure  6:  Example  of  interpolation  weighting  factors  based  on  north-south 
distance.  In  this  example,  the  interpolation  is  from  NAVDAS  grid  points  to  the 
NOGAPS  grid  point  at  27.7N,  which  represents  a  distance  of  1.5°  from  26.95N  to 
28.45N.  The  total  weighting  factor  depends  also  on  the  east-west  grid 
configuration  (Fig.  7). 
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Weighting  factors  for  horizontal  interpolation  (like  those  shown  in  Fig.  5)  are 
determined  by  a  linear  algorithm  based  on  distances  in  the  north-south  and  east- 
west  directions.  Consider  the  ’’forward”  interpolation  from  the  NAVDAS  grid  to 
the  NOGAPS  grid.  Each  T79  grid  point  represents  a  north-south  distance  of  1.5° 
(approximately  1.5°  on  the  actual  T79  grid),  centered  on  the  gridpoint  latitude.  If 
the  latitude  of  the  T79  gridpoint  is,  say,  27.7N,  it  represents  the  interval  from  26.95N 
to  28.45N,  as  shown  in  Fig.  6.  This  interval  on  the  T79  grid  is  overlapped  by  intervals 
represented  by  several  (either  two  or  three)  points  on  the  NAVDAS  1°  grid.  In  the 
example  of  Fig.  6  the  interpolation  to  the  T79  grid  point  at  27.7N  is  a  weighted 
average  of  values  taken  from  points  on  the  NAVDAS  1°  grid  at  27N  and  28N.  The 
total  weighting  factor  for  horizontal  interpolation  is  a  product  of  the  north-south 
weighting  factor  and  another  factor  that  depends  on  the  east- west  grid  configuration, 
as  shown  in  Fig.  7.  If  the  longitude  of  the  T79  gridpoint  is  an  integer  (e.g.,  it  lies 
directly  on  the  1°  grid  at  gridpoint  ”i”)  the  assigned  weighting  is  1/6  to  the  column 
indexed  (i-1),  2/3  to  the  column  indexed  (i)  and  1/6  to  the  colmnn  indexed  (i+1). 
In  the  alternate  case  where  the  T79  gridpoint  is  halfway  between  points  on  the  1° 
grid,  the  two  columns  to  the  east  and  west  are  each  assigned  half  the  weight. 
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Longitude 


■  NOGAPS  T79  grid 


G  NAVDAS  rgrid 


Figure  7:  Interpolation  weighting  based  on  longitudinal  position  of  NAVDAS  and 
NOGAPS  gridpoints.  The  total  horizontal  interpolation  weighting  factor  as  shown 
for  selected  gridpoints  in  Fig.  5  combines  a  weighting  factor  for  the  east-west 
direction  (illustrated  here)  and  a  weighting  factor  for  the  north-south  direction  as 
illiKtrated  in  Fig.  6. 


3.2  Output  File 


The  output  of  the  NAVDAS  adjoint  is  written  to  a  formatted  file.  Each  Une  of  the  out¬ 
put  file  corresponds  to  a  separate  observation  (excluding  any  rejected  observations), 
and  includes  the  following  information: 
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•  n  -  sequence  number  of  observation 


•  rlat_ob  -  latitude  of  observation 

•  rlon_ob  -  longitude  of  observation 

•  ob_sens  -  observation  sensitivity,  dJ  /  dy 

•  adjoint  cob  -  (  z  =  H  Pj,  a  J  /  dxa) 

•  jvarty_ob  -  observation  type  (l=height,  2=temperature,  3=u-wind,  4=v-wind, 
13=brightness  temperature) 

•  insty_ob  -  instrument  type 

•  num  reject  -  reject  flag  (0=assimilated,  l=rejected) 

•  err_ob  -  assigned  observation  error  standard  deviation 

•  rmsvar  -  assigned  background  error  standard  deviation 

•  P_ob  -  pressure  level  of  observation 

•  c_pf_ob  -  observation  label  (pt.  1) 

•  c_db_ob  -  observation  label  (pt.  2) 

•  xiv_ob  -  innovation  (observation  -  background) 

•  ob  -  observation  value 


26 


•  oberr_sens  -  sensitivity  to  observation  error  standard  deviation 

•  rsd_val  -  residual  (analysis  -  observation) 


3.3  Gradient  (accuracy)  Test 


A  basic  procedure  for  validating  that  an  adjoint  model  has  been  correctly  coded  is  the 
so-called  gradient  test,  which  compares  the  value  of  scalar  inner  products,  using  input 
and  output  sensitivity  vectors.  With  the  NAVDAS  adjoint,  a  form  of  the  gradient 
test  can  be  defined  as 


(y  -  Hxfo)  , 


dy)  V  '  ’ 


(14) 


observation  —  space 
^Johs 


gridpoint  —  space 
dJgrid 


where  the  LHS  of  (14)  is  calculated  in  observation  space  using  temperature,  wind,  and 
height  innovations  and  sensitivity,  and  the  RHS  is  calculated  in  the  gridpoint  space  of 
the  forecast  model  using  temperature,  wind,  and  terrain  pressure  anal3^is  increments 
and  sensitivity.  For  consistency  with  results  shown  later,  it  is  convenient  to  define 
the  costfunction  here  so  that  5J  represents  the  difference  between  the  forecast  errors 
of  trajectories  starting  from  Xo  and  x^.  Thus,  e/  is  the  error  in  a  forecast  of  length  / 
starting  fi:om  Xa,and  Cg  is  the  error  in  a  forecast  of  length  g  starting  from  x;,.  The  6-hr 
forecast  of  trajectory  g  provides  the  background  for  the  assimilation  that  produces 
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Xa,  and  the  two  forecasts  verify  at  the  same  time.  The  scalar  quantities  5  Jobs 
SJgrid  ^re  adjoint-based  estimates  of  the  actual  error  difference  Bf  —  eg.  Appendix  B 
provides  a  derivation  of  the  expressions  used  in  (14).  Note  that  the  forecast  length 
does  not  affect  the  accuracy  of  this  gradient  test,  since  the  test  only  depends  on  the 
accuracy  of  the  NAVDAS  adjoint  calculation  (e.g.,  Eq.  1)  and  not  on  the  accuracy 
of  the  forecast  model  (NOGAPS)  adjoint.  However,  forecast  length  does  affect  the 
accuracy  of  5  Jobs  and  dJgru  with  respect  to  the  actual  value  of  e/- eg  that  is  obtained 
from  the  nonlinear  forecast  model  (see  Section  4.2). 


Figure  8:  Gradient  test  comparison  for  &/obs  (solid  line)  and  &/grid  (dash  line)  for 
Dec  2002.  Costfunction  J  is  NOGAPS  forecast  error  (624-63^  ,  J  kg'')  in  the 
global  domain.  St/obs  and  &/gridare  defined  in  Section  3.3. 
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A  comparison  of  SJobs  with  dJgrid  using  /  =  2Ah  and  g  =  ZOh  for  the  month  of 
December  2002  is  shown  in  Fig.  8.  The  forecast  errors  are  calculated  in  the  global 
domain  from  NOGAPS  forecasts  started  at  OOUTC.  For  December  2002,  the  average 
value  of  5Jobs/SJgrid  is  0.951  (maximum  =  1.002,  minimum  =  0.883).  These  results 
indicate  that  the  NAVDAS  adjoint  calculates  the  observation  sensitivity  from  the 
analysis  sensitivity  input  to  within  about  95%  accuracy.  Note  that  6J  is  negative, 
since  624  <  630. 

In  theory  the  equivalence  of  SJobs  and  SJgrid  should  be  exact  if  the  adjoint  code 
that  produces  the  sensitivity  dJ  /  dy  is  completely  consistent  with  respect  to  the 
version  of  NAVDAS  that  produces  the  analysis  increment  (xq— xt)  from  the  innovation 
(y — Hxfc).  In  practice,  there  are  several  obstacles  that  prevent  an  exact  result.  First, 
the  interpolation  of  dJ  /  5xa  from  NOGAPS  to  NAVDAS  gridpoint  space  unavoidably 
introduces  some  small  (but  finite)  inaccuracy  into  the  observation  sensitivity  results. 

An  additional  consideration  is  that  the  NOGAPS  adjoint  produces  a  sensitivity 
to  the  analyzed  terrain  pressure  as  a  component  of  dJ  /  5xo,  but  this  pressure 
sensitivity  does  not  directly  translate  into  a  component  of  the  sensitivity  input  to 
NAVDAS,  because  NAVDAS  does  not  produce  a  pressine  analysis  increment.  Since 
the  terrain  pressure  contribution  to  the  input  analysis  sensitivity  from  NOGAPS  is 
neglected  in  Fig.  8,  there  is  a  small  systematic  underestimate  of  the  observation 
sensitivity  magnitude  in  the  NAVDAS  adjoint  calculation. 
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3.4  Code  Maintenance  and  Updating 


The  current  version  of  the  NAVDAS  adjoint  is  consistent  with  the  forward  version 
of  NAVDAS  as  of  January  2004.  It  is  anticipated  that  the  adjoint  code  will  be 
updated  at  periodic  intervals  in  order  to  maintain  applicabihty  to  the  operational 
data  assimilation  system.  For  example,  new  observation  quantities  such  as  ozone  will 
be  added  to  the  adjoint  code.  However,  the  NAVDAS  adjoint  should  be  considered 
a  research  code  and  it  is  not  part  of  the  formal  NAVDAS  configuration  management 
system.  Refer  to  Appendix  A  for  the  locations  of  current  codes  and  scripts. 

4  Observation  Sensitivity  Examples 

The  current  version  of  the  NAVDAS  adjoint  provides  two  sensitivity  gradients,  i)  the 
sensitivity  to  observations  dJ  /  dy  and  ii)  the  sensitivity  to  the  assigned  observation 
error  (see  Section  4.5).  The  observation  sensitivity  dJ  f  dy  can  be  used  for  a  variety 
of  applications,  which  involve  various  choices  of  observation  (or  innovation)  pertur¬ 
bation  6y  and  costfunction  J.  Four  possible  applications  are  summarized  in  Table  2. 
Application  ^1  is  the  use  of  observation  sensitivity  to  interpret  arbitrary  perturba¬ 
tions  of  the  observations  or  background.  In  application  ^2,  we  quantify  the  impact 
of  actual  innovations  on  a  known  forecast  error  difference.  Application  #3  describes 
the  general  interpretation  of  sensitivity  for  hypothetical  observations.  Application 
^4  pertains  to  targeted  observing,  when  the  both  innovations  and  forecast  errors  are 
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not  known. 


Application 

8y 

/(scalar) 

1.  Basic  Observation  Sensitivity: 

What  are  the  impacts  of  various 
observation  perturbations  ? 

Any  perturbation 
of  appropriate 
size 

Any  valid 
costfunction 

2.  Innovation  Impact: 

What  are  the  impacts  of  actual  innovations 
on  a  known  forecast  error  ? 

y  -  Hxb  (all 
knovra) 

Forecast  error 
(known) 

3.  Hypothetical  Observations: 

What  are  the  impacts  of  hypothetical 
observations  on  a  known  forecast  error  ? 

y  -  Hxb  (some 
known,  some 
imknown) 

Forecast  error  or 
other  costfunction 

4.  Targeted  Observing: 

What  is  the  impact  of  adding  a  set  of 
targeted  observations  to  the  regular 
observing  network  at  a  future  time? 

y  -  Hcb  (all 
unknown) 

Costfunction 
usually  a  forecast 
error  surrogate 

Table  2:  Examples  of  observation  sensitivity  applications. 
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4.1  Basic  Observation  Sensitivity 


In  the  most  general  application,  observation  sensitivity  can  be  used  to  estimate  the 
effect  of  various  arbitrary  changes  in  the  values  of  observations  (or  innovations).  For 
example,  we  can  define  J  as  the  difference  between  the  24h  and  30h  global  forecast 
error  and  examine  the  sensitivity  of  this  costfunction  with  respect  to  any  component 
of  the  observation  vector,  as  illustrated  in  Fig.  9  for  rawinsonde  temperature  obser¬ 
vations  assimilated  at  OOUTC  10  December  2002.  The  sensitivity  vector  dJ  /  dy  is 
plotted  in  observation  space  at  the  locations  where  observations  have  been  assimi¬ 
lated  by  NAVDAS.  We  can  also  define  an  observation  perturbation  vector  ^y,  which 
can  represent  a  change  to  the  value  of  any  one,  or  all,  observations  assimilated  by 
NAVDAS.  Using  (12),  the  perturbation  6y  implies  that  the  forecast  costfunction  will 
be  changed  by  an  amount  SJ.  Thus,  if  we  consider  one  temperature  observation  for 
which  (suppose)  dJ  j  dy  ^  0.004  J  kg-Meg-^  and  we  let  by  ^  2.0  deg,  then  bJ 
—  0.008  J  kg  The  magnitude  of  by  is  normally  restricted  to  perturbations  no 
larger  than  typical  innovations  because  of  tangent  linear  approximations  in  the  fore¬ 
cast  model  adjoint,  and  the  sensitivity  information  is  most  accurate  for  short-range 
forecasts  of  72h  or  less. 

In  addition  to  estimation  of  observation  perturbation  impact,  the  observation 
sensitivity  dJ  f  dy  and  the  analysis  sensitivity  dJ  /  dxa  can  be  used  to  determine 
the  value  of  bJ  implied  by  a  background  perturbation  using  (13),  as  described 
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in  Section  2.2. 


Figure  9:  Sensitivity  of  J  =  e24-e3(,  to  500  mb  rawinsonde  temperature  observations  at 
OOUTC  10  Dec  2002.  Units  are  lO-^  J  kg->deg-'. 


4.2  Observation  Impact  On  Forecast  Error 

Observation  sensitivity  can  be  used  in  a  diagnostic  mode  to  evaluate  the  impact 
of  actual  observations  (or  innovations)  on  the  difference  between  two  short-range 
forecast  errors.  As  in  Section  3.3,  we  define  e/  as  the  error  in  a  forecast  of  length 
/  starting  from  Xa,  and  eg  is  the  error  in  a  forecast  of  length  g  starting  from  x^. 
The  6-hr  forecast  of  trajectory  g  provides  the  background  for  the  assimilation  that 


33 


produces  x^,  and  the  two  forecasts  verify  at  the  same  time.  In  general,  Xq  will  be 
a  better  estimate  of  the  true  atmospheric  state  than  Xb,  and  thus  typically  (but  not 
always)  624  <  630-  The  error  difference  e/  —  eg  is  due  entirely  to  the  assimilation  of 
observations  used  to  produce  the  analysis  Xq.  The  true  difference  between  the  forecast 
errors  e/  and  eg  is 

Ae®  =  ef-eg.  (15) 

Note  that  in  a  limiting  case  in  which  no  observations  are  assimilated,  Xq  =  xt  and 
Ae®  =0.  As  shown  in  Langland  and  Baker  (2004)  and  in  Appendix  B,  observation 
sensitivity  can  be  used  to  provide  an  estimate  of  Ae^,  using  what  can  be  called  an 
’’observation  impact”  equation 

Se}^({y-Hxb)  ,  .  (16) 

For  the  global  domain,  the  inner  product  represented  by  (16)  is  evaluated  using  the 
entire  innovation  and  observation  sensitivity  vectors  (excluding  moisture  at  present). 
Note  that  for  this  choice  of  costfunction  Se^j  in  (16)  is  the  same  quantity  as  5  Jobs 
in  (14).  It  can  also  be  noted  that  studies  by  Doerenbecher  and  Bergot  (2001)  and 
Fourrie  et  al.  (2002)  use  observation  impact  functions  in  observation  space  that 
are  similar  in  form  to  Elq.  (16)  -  e.g.,  they  involve  the  innovation  vector  and  an 
observation  sensitivity  gradient  However,  in  those  studies  the  observation  sensitivity 
is  not  derived  with  an  actual  adjoint  of  the  assimilation  procedure.  Fourrie  et  al. 
(2002)  use  the  error  difference  e/  —  as  a  costfunction,  but  their  impact  function  is 
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based  on  sensitivity  gradients  that  involve  only  the  trajectory  starting  from  Xq.  In 
Doerenbecher  and  Bergot  (2001),  the  costfunction  is  a  dynamic  variable  (enstrophy), 
rather  than  a  forecast  error. 


Figure  10:  Time  series  of  624-630  (dark  solid  line)  in  NOGAPS 
forecasts  and  corresponding  adjoint-based  observation  impact 
estimate  (grey  solid  line)  calculated  using  Eq.  16  in  the  global 
domain  with  no  moisture  observations  in  December  2002. 


A  time  series  comparing  Ae®  and  Se^  for  December  2002,  using  /  =  24h  and 
g  =  30/i  is  provided  in  Fig.  10.  The  values  of  5e^  and  Ae^®  are  negative,  since 
624  <  630.  For  December  2002,  the  average  value  of  5e^  /  Aej®  is  0.740  (maximum 
=  0.827,  minimum  =  0.672).  The  observation  impact  calculations  using  (16)  are 
thus  an  underestimate  of  the  actual  forecast  error  difference.  The  majority  of  the 
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underestimate  is  probably  explained  by  the  neglect  of  moisture  observations^  in  the 
current  NAVDAS  adjoint.  However,  the  accuracy  of  the  observation  impact  estimate 
is  also  affected  by  approximations  in  the  forecast  model  adjoint,  primarily  the  neglect 
of  moist  physical  processes  and  nonlinearities  in  the  dry  dynamics.  The  effect  of  these 
tangent  linear  approximations  becomes  more  significant  as  forecast  length  increases. 


Ob  count 

3,000,000 

2,000,000. 

1.000,000 


Figure  11:  Summed  global  observation  impact  (  6e^ ,  J  kg-^)  for 
Southern  and  Northern  Hemisphere,  partitioned  by  instrument 
type,  combining  June  and  December  2002.  Includes  all 
observations  assimilated  in  NAVDAS  at  OOUTC.  ATOVS- 
temperature  relrievals,  RAOB-rawinsondes,  SATW-cloud  and 
feature-track  winds,  AIRW-commercial  aircraft  observations, 
LAND-land  surface  observations,  SHIP-ship  surface  observations, 
AUSN- synthetic  sea-level  pressure  data 


^Moisture  data  account  for  about  30  percent  of  total  global  observations  at  OOUTC. 


36 


Fig.  11  illustrates  how  observation  impact  estimates  can  be  used  to  compare  the 
value  provided  by  the  various  types  of  observations  assimilated  by  NAVDAS.  Here, 
observation  impact  is  evaluated  in  terms  of  Sel^  for  the  northern  and  southern  hemi¬ 
spheres,  combining  the  months  of  June  and  December  2002  and  using  observations 
assimilated  at  OOUTC.  ATOVS  are  satellite  temperature  retrievals  (Reale  et  al.  2003). 
SATW  are  feature-track  wind  vectors  from  visible,  microwave  and  water  vapor  geo¬ 
stationary  imagery  (Rao  et  al.  2002).  RAOB  are  rawinsonde  temperatmre,  winds,  and 
surface  heights.  AIRW  are  commercial  aircraft  temperature  and  wind  data  in  level 
flight  and  ascent  and  descent  profiles.  LAND  are  surface  temperatures,  winds,  and 
heights  from  land  surface  reporting  stations.  SHIP  are  surface  temperatures,  winds, 
and  heights  from  ships  at-sea.  AUSN  are  synthetic  sea-level  pressure  data  assimilated 
as  height  observations  over  part  of  the  southern  hemisphere  (Guymer  1978). 

As  shown  in  Fig.  11,  the  largest  forecast  error  reductions  (Sell)  for  the  south¬ 
ern  hemisphere  from  observations  assimilated  at  OOUTC  are  produced  by  ATOVS, 
satellite  wind  data,  and  rawinsondes.  In  the  northern  hemisphere  the  largest  error 
reductions  are  produced  by  rawinsonde,  satellite  wind  data,  commercial  aircraft  data 
and  ATOVS.  Almost  all  the  error  reduction  due  to  commercial  aircraft  observations 
is  attributable  to  data  in  the  northern  hemisphere.  Land  and  ship  siuface  data  pro¬ 
duce  lesser  (but  still  significant)  forecast  error  reductions.  AUSN  data  are  a  valuable 
observation  type  for  the  southern  hemisphere.  In  fact,  the  error  reduction  per  obser- 
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vation  for  AUSN  can  be  much  larger  than  for  other  types  of  observations  (Langland 
and  Baker  2004). 


a  Rawinsonde  Profiles  b  ATOVS  Profiles 


Figure  12:  Observation  impact  (J  kg-')  in  the  NAVDAS  assimilation 
at  OOUTC  10  December  2002.  Green  (red)  dots  represent  large 
reduction  (increase)  in  24h  global  forecast  error.  Blue  (orange)  dots 
represent  moderate  reduction  (increase)  in  24h  global  forecast  error. 
Grey  dots  represent  relatively  small  reduction  or  increase  in  24h  global 
forecast  error.  For  rawinsondes  and  ATOVS  each  dot  represents  the 
combined  impact  of  observations  in  vertical  profiles.  For  aircraft  and 
satellite  wind  observations  each  dot  combines  the  impact  of 
observations  with  the  same  latitude,  longitude  and  pressure  level.  Z  = 
global  sum  of  observation  impact  for  instrument  type. 
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While  Fig.  11  indicates  that  the  global  (or  hemispheric)  sums  of  5e^  are  negative 
(confirming  that  assimilation  of  large  numbers  of  observations  reduces  624  —  630)  these 
sums  combine  large  numbers  of  positive,  as  well  as  negative,  observation  impacts 
associated  with  individual  observations.  This  can  be  seen  by  examining  the  pattern  of 
observation  impact  at  a  particular  assimilation  time.  For  example,  at  the  assimilation 
time  of  OOUTC  10  Dec  2002  (Fig.  12),  each  of  the  four  observation  types  include 
many  individual  data  for  which  Se^  >  0  (shown  as  orange  and  red  dots).  This  mix 
of  positive  and  negative  observation  impact  at  a  single  assimilation  time  does  not 
indicate  a  systematic  problem  with  certain  observations,  rather  it  is  an  outcome  of 
the  statistical  assumptions  and  approximations  used  in  data  assimilation.  However, 
if  the  impact  of  a  particular  observation  type  or  reporting  station  is  consistently  > 
0,  it  could  indicate  an  instrmnent  or  quality  control  problem. 
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Top  20  Raobs  -  24h  Post  Error 


- — ^ 
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Figure  13:  Twenty  rawinsonde  observing  stations  that  produce  the 
largest  observation  impact  Se^  ,  representing  reductions  of  624-030 
the  global  domain,  due  to  observations  assimilated  at  OOUTC  in 
December  2002. 


Another  example  of  observation  impact  (Fig.  13)  shows  the  locations  of  20  rawin¬ 
sonde  stations  that  produced  the  largest  negative  values  (e.g.,  error  reduction)  of 
for  OOUTC  assimilations  and  forecasts  during  the  month  of  December  2002.  These 
key  rawinsondes  tend  to  be  found  along  land  /  ocean  boundaries,  in  regions  where 
forecast  error  growth  rates  are  relatively  large,  and  where  other  in-situ  or  satellite 
observations  are  relatively  sparse. 
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4.3  Hypotheticed  Observations 


The  NAVDAS  adjoint  can  be  used  to  determine  the  sensitivity  of  a  known  forecast 
error  or  forecast  parameter  to  hypothetical,  as  well  as  "real”  observations.  The  fol¬ 
lowing  information  must  be  provided  for  each  of  the  hypothetical  observations: 

•  Latitude  and  Longitude 

•  Pressure  Level 

•  Instrument  Type  (e.g.,  rawinsonde,  satellite  wind,  etc.) 

•  Observation  type  (e.g.,  temperature,  wind,  height) 

•  Specified  observation  error 

In  this  application  the  costfunction  ( J)  can  be  forecast  error  or  some  other  fore¬ 
cast  parameter  of  interest,  such  as  surface  pressure,  or  wind  speed.  The  NAVDAS 
adjoint  will  produce  sensitivity  to  both  the  real  and  hypothetical  observations.  Note 
that  the  observation  sensitivity  can  be  determined  without  the  specification  of  an 
innovation  value.  The  impact  of  the  real  observations  on  an  actual  short-range  fore¬ 
cast  error  difference  can  be  calculated  using  (16).  The  impact  of  the  hypothetical 
observations  can  also  be  estimated  using  (16),  if  a  proxy  innovation  value  is  specified. 
Sensitivity  to  hypothetical  observations  can  be  used,  for  example,  to  develop  more 
optimal  configurations  of  satellite  and  in-situ  observing  systems. 
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Figure  14:  Configuration  for  targeted  observing  used  in  Fig.  15. 
The  interval  from  +0h  to  +66h  is  used  to  produce  the  targeting 
guidance  and  deploy  observational  resources  to  the  target  area. 


4.4  Tao'geted  Observing 

Observation  sensitivity  can  be  used  for  ’’targeted  observing,”  in  which  objective 
(model-based)  information  provides  guidance  for  the  deployment  of  special  obser¬ 
vations  intended  to  improve  a  short-range  numerical  weather  forecast  over  a  specific 
region.  This  is  a  variant  of  sensitivity  to  hypothetical  observations,  in  which  a 
costfunction  other  than  actual  forecast  error  must  be  used,  and  the  observation  sen¬ 
sitivity  is  calculated  for  a  future  time  (typically  24-72h  ahead)  when  the  targeted 
observations  are  to  be  assimilated.  An  example  of  a  possible  time  configuration  for 
targeted  observing  is  shown  in  Fig.  14.  Usually,  the  ’’targeted  observations”  will  be 
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a  relatively  small  set  of  special  data  added  to  the  observations  provided  by  regular 
satellite  and  in-situ  observing  systems.  The  goal  of  targeted  observing  is  to  maximize 
the  value  of  the  complete  observing  network,  e.g,  the  regular  observations  plus  the 
special  targeted  observations.  However,  adding  targeted  observations  will,  in  general, 
change  the  impact  of  nearby  regular  observations.  The  targeted  observing  procedure 
generally  includes  the  following  steps: 

Procedure  for  Targeted  Observing: 

1.  Identify  the  forecast  feature  that  is  to  be  improved,  through  examination  of 
deterministic  or  ensemble  forecasts 

2.  Define  a  forecast  verification  region,  verification  time  and  costfunction  ( J) 

3.  Define  the  type  of  targeted  observation(s)  and  time  at  which  they  are  to  be 
obtained 

4.  Define  several  (or  many)  possible  configurations  for  the  targeted  observations 

5.  Calculate  observation  sensitivity  dJ  /  dy  for  each  configuration  (each  requires 
a  separate  run  of  NAVDAS  adjoint) 

6.  Prioritize  targeting  configiuration  based  on  impact  fimction  (e.g.,  Eq.  18,  Eq. 
19,  or  other  expression) 

7.  Obtain  targeted  observations  using  in-situ  or  satellite  instruments 

8.  Assimilate  targeted  observations  and  produce  forecast 
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In  step  5,  the  locations  of  the  hypothetical  targeted  observations  can  be  spec¬ 
ified  using  an  input  file  read  by  the  NAVDAS  adjoint.  The  locations  of  regular 
rawinsonde,  land-surface,  and  ATOVS  observations  at  the  future  targeting  time  can 
be  estimated  with  fairly  good  accuracy.  The  locations  of  commercial  aircraft  and 
feature-track  satellite  wind  data  can  be  approximated  using  information  from  the 
most-recent  analysis  performed  at  the  same  hour  (e.g.,  OOUTC,  18UTC,  etc.)  as  the 
tune  at  which  the  targeted  observations  are  to  be  taken.  The  costfunction  J  for  tar¬ 
geted  observing  can  be  a  forecast  parameter,  such  as  vorticity,  wind,  temperature,  or 
surface  pressure,  or  it  can  be  a  difference  between  two  forecasts,  but  the  costfunction 
in  targeted  observing  cannot  be  forecast  error,  since  that  is  not  known  in  advance. 
Another  constraint  is  that  the  actual  innovation  values  are  not  known  for  any  of  the 
observations,  either  real  or  hypothetical,  that  will  be  obtained  at  the  targeting  time. 
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Figure  15:  Estimated  impact  of  two  dropsonde  configurations  on  the 
42h  NOGAPS  forecast  for  northern  Europe  (area  within  red  dash  line). 
Targeted  observing  time  is  18UTC  15  Nov  2003,  forecast  verification 
time  is  12UTC  17Nov2003. 


Fig.  15  illustrates  the  results  of  an  observation  sensitivity  calculation  for  two 
possible  targeting  observing  configurations  involving  dropsondes  deployed  by  a  re¬ 
connaissance  aircraft.  In  this  case,  the  targeted  observing  time  is  18UTC  15  Dec 
2003  and  the  forecast  verification  time  is  12UTC  17  Dec  2003  (-b42h).  The  forecast 
verification  region  is  an  area  of  the  North  Atlantic  and  northern  Europe  (indicated 
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by  the  dashed  red  box).  The  NOGAPS  forecast  trajectory  used  for  the  targeting 
calculations  starts  from  an  analysis  at  OOUTC  13  Dec  2003.  Thus,  in  this  example, 
there  is  a  66h  ”lead  time”  in  which  to  produce  the  targeting  guidance  and  deploy 
an  aircraft  to  the  target  region.  The  cost  function  (a  proxy  for  the  actual  forecast 
error  or  forecast  uncertainty)  used  here  is  the  energy- weighted  difference  between  two 
forecasts  that  verify  at  12UTC  17  Dec  2003  -  a  108h  forecast  from  OOUTC  13  Dec 
2003  and  a  114h  forecast  from  18UTC  12  Dec  2003. 

^  2  ^  (^108  “  X114)  ,  C  (Xi08  —  X114))  ,  (17) 

where  C  is  a  matrix  of  energy  weighting  coefficients  for  dry  total  energy  (Rosmond 
1997),  and  the  NOGAPS  forecast  state  vector  x  includes  temperature,  vorticity,  diver¬ 
gence,  and  surface  pressure.  Note  that  for  targeted  observing,  it  is  not  known  whether 
the  desired  outcome  (a  reduction  in  forecast  error)  corresponds  to  an  increase  or  de¬ 
crease  in  the  value  of  the  costfunction.  An  improved  forecast  might  correspond  to 
either  an  increase  or  decrease  in  forecast  energy,  or  surface  pressure  over  the  verifica¬ 
tion  region,  for  example.  Thus,  the  impact  function  for  targeted  observing  can  only 
measure  the  expected  unpact  of  the  observations  in  terms  of  a  change  in  the  variance 
or  magnitude  of  J,  under  the  assumption  that  if  the  observations  have  an  impact 
on  J  there  is  also  a  potential  for  reduction  in  forecast  error.  One  expression  for  the 
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expected  impact  of  a  configuration  (n)  of  targeted  observations  is 

regular  _ohs  taTgeted_ohs 

l<^‘^t{n)  =  (  |y  -  Hxbl  .  )  +  (  |y  -  Hxa!  .  )  •  (18) 

\^sumed  mnovation  j  \assumed  mnovation  j 

Eq.  (18)  is  a  generalization  of  the  observation  impact  equation  (16)  for  the  sit¬ 
uation  when  the  iimovations  and  actual  forecast  errors  are  not  known.  Since  the 
innovations  are  not  known  when  the  targeting  guidance  is  prepared,  it  is  necessary 
to  use  ’’assumed  innovations”  in  (18).  The  assumed  innovation  can  be,  for  example, 
an  estimate  of  analysis  uncertainty  or  background  error,  or  even  a  constant,  in  which 
case  is  directly  proportional  to  the  sensitivity  gradient  magnitude.  In  (18) 

both  the  ’’innovation”  and  the  observation  sensitivity  appear  as  absolute  values,  so 
that  each  observation  contributes  to  a  positive  value  of  .  The  total  value  of 

is  the  inner  product  calculated  over  the  global  domain,  including  all  regular 
and  targeted  observations. 

In  Fig  15a,  Flight  #1  is  a  set  of  20  h3q)othetical  dropsondes  deployed  from  an 
altitude  of  24,000  ft.  The  location  of  each  dropsonde  is  indicated  by  a  solid  dot, 
and  each  dropsonde  provides  a  vertical  profile  of  temperature  and  wind  on  pressure 
surfaces  at  intervals  of  100  mb,  with  the  assumed  accuracy  of  rawinsonde  observations. 
Flight  #2,  shown  in  Fig.  15b,  is  a  set  of  12  hypothetical  dropsondes  with  the  same 
specifications  as  Flight  #1.  The  other  solid  dots  in  Fig.  15  are  locations  of  regular 
rawinsonde  stations  providing  observations  at  the  18UTC  time. 
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Observation 

..Type _ 

Control 

Flight  #1 

Flight  #2 

Raobs 

0.1520 

0.1475 

0.1483 

Dropsondes 

- 

0.0461 

0.0578 

Satwinds 

0.1560 

0.1530 

0.1536 

Aircraft 

0.3839 

0.3794 

0.3832 

ATOVS 

0.5003 

0.4911 

0.4894 

Land  Surface 

0.0144 

0.0143 

0.0143 

Ship  Surface 

0.0695 

0.0690 

0.0668 

Total 

1.2761 

1.3004 

1.3134 

Table  3:  Impact  function  |&/|  (J  kg’')  for  targeted  observing  flight  tracks 
shown  in  Fig.  15.  Costfimction  (J)  is  the  energy-weighted  difference  between 
108h  and  114h  forecasts  verifying  over  Northern  Europe  at  12UTC  17  Dec 
2003. 


Table  3  summarizes  the  targeting  impact  calculations  using  (18)  with  an  assumed 
innovation  of  1.0  and  observation  sensitivity  calculated  with  the  NAVDAS  adjoint. 
With  no  dropsondes  added,  the  value  of  |5J|  for  this  case  is  equal  to  1.2761  J  kg“^ 
This  number  indicates  the  ’’potential”  of  the  regular  observing  system  to  influence 
the  forecast  uncertainty  represented  by  the  costfimction.  If  no  observations  are  as¬ 
similated,  then  |5J|  =  0. 

By  adding  various  configurations  of  targeted  observations,  we  can  increase  the 
potential  impact  of  the  total  observing  system.  The  20  dropsondes  in  Flight  #1 
increase  |5i7|  to  1.3004  J  kg  ^  (Table  3).  Note  that  adding  the  dropsondes  decreases 
marginally  the  impact  of  other  observation  types,  but  increases  the  potential  impact  of 


48 


the  total  (regular  +  targeted)  observing  network,  which  is  the  desired  result.  Flight 
#2  includes  only  12  dropsondes,  but  increases  |(5J|  to  1.3134  J  kg“^  and  is  thus 
considered  to  be  the  better  of  the  two  possible  dropsonde  configurations.  In  fact, 
Flight  t^2  covers  an  area  of  stronger  analysis  sensitivity  (not  shown).  This  targeting 
example  suggests  a  relatively  modest  forecast  impact,  changing  |5J|  by  approximately 
2-3%.  However,  during  the  Winter  Storm  Reconnaissance  Program  (WSR)  targeted 
dropsonde  observations  reduced  short-range  forecast  errors  by  an  average  10%-20% 
(Szimyogh  et  al.  2000)  and  by  up  to  50%  over  localized  regions  during  the  North 
Pacific  Experiment  (NORPEX,  Langland  et  al.  1999). 

An  alternative  equation  to  estimate  the  impact  of  targeted  observations  can  be 
obtained  (following  Baker  2000)  if  we  define  J  as  the  projection  of  the  analysis  error 
Sa  onto  the  analysis  sensitivity  gradient 

.  (w) 

The  expected  variance  of  the  change  in  the  forecast  aspect  J  is 

where  )  is  the  analysis  error  covariance  matrix,  given  by 

Pa  =  Pb-  (HPfeH^  -h  R)  HPfe  .  (21) 

The  second  term  in  (21)  represents  the  reduction  of  the  backgroimd  error  covari¬ 
ance  due  to  the  presence  of  the  observations.  Eqs  (20)  and  (21)  can  be  combined 
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(22) 


(Baker  2000),  finally  resulting  in 

((")’>=  (|)"|HP.H-  +  R|  I . 

In  (22)  {{SJf)  is  the  expected  reduction  in  the  variance  of  5  J  due  to  the  targeted 
and  regular  observations.  This  expression  can  be  evaluated  using  the  observation 
sensitivity  and  ”cob”  vectors,  both  of  which  are  provided  by  the  NAVDAS  adjoint. 

4,5  Sensitivity  to  Observation  Error 

The  quality  of  a  data  assimilation  procedure  depends  partly  on  the  specification 
of  observation  and  background  error.  These  error  parameters  are  generally  tuned  or 
adjusted  based  on  information  from  various  sources,  which  may  include  values  used  in 
other  data  assimilation  procedures  or  the  results  of  forward  ’’sensitivity”  experiments 
in  which  specified  errors  are  modified  in  more-or-less  arbitrary  ways.  The  trial  and 
error  sensitivity  approach  is  generally  an  inefficient  tuning  procedure  and  is  unlikely 
to  produce  the  optimal  specification  of  either  observation  or  background  error.  For 
example,  the  results  of  conventional  sensitivity  tests  are  often  ambiguous  because  the 
’’tuning”  generally  produces  a  combination  of  positive  and  negative  impacts  that  are 
partially  self-canceling  and  conceal  the  actual  potential  for  larger  impact. 

A  potentially  more  efficient  method  of  parameter  tuning  can  use  adjoint  sensitivity 
information  which,  for  a  selected  forecast  costfunction,  provides  a  complete  sensitivity 
gradient  vector  for  the  entire  observation  space  of  the  data  assimilation  procedure. 
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With  this  information,  the  specified  observation  errors  (or  other  parameters)  can  be 
adjusted  to  achieve  more  optimal  results  with  less  trial  and  error. 

The  values  of  the  observation  error  standard  deviation  £j.  in  NAVDAS  are  specified 
by  data  statements  in  the  innovation  code  module  and  written  to  the  innovation 
file.  In  the  NAVDAS  adjoint  code  the  innovation  file  is  read  and  elements  of  the 
diagonal  observation  error  covariance  matrix  R  are  defined  as  ri  =  (Sr  /  where 
Sb  is  the  background  error  standard  deviation  at  observation  locations,  and  Sb  is  used 
here  as  a  normalization  factor.  We  wish  to  obtain  the  sensitivity  gradient  dJ  /  dvi, 
which  will  be  a  vector  in  observation  space.  One  method  to  obtain  dJ  /  dvi  would 
be  to  write  an  adjoint  (transpose)  version  of  the  discretized  code  that  defines  the 
NAVDAS  operator  +  R]“^  in  which  £r  and  Ti  appear.  This  has  not  been 

done  because  the  operator  is  self-adjoint  and  additional  adjoint  code  development 
was  not  necessary  to  obtain  sensitivity  to  observations. 

Alternatively,  it  is  possible  to  define  dJ  /  dri  (the  array  oberr_sens)  using  vari¬ 
ables  already  provided  by  the  forward  and  adjoint  versions  of  NAVDAS  as 

dJ  dJ 
dn~ 

where  z  is  the  ”cob”  array  defined  by  the  solver  in  the  forward  NAVDAS  ( z  =[HP6H'^-|- 
R]  (y  —  Hxft)).  The  NAVDAS  adjoint  reads  z  from  the  rsd_vec  file  and  then  com¬ 
putes  dJ  j  dvi  after  the  calculation  oidJ  /  dy  has  been  completed.  The  s3unbol  (*) 
in  (20)  denotes  a  vector  Shur  product,  not  a  vector  inner  or  outer  product.  Note  that 
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dJ  /  dri  has  the  units  of  J,  since  ri  is  non-dimensional.  The  array  z  obtained  from 
the  rsd_vec  file  has  the  units  of  y.  The  derivation  of  (20)  appears  in  Appendix  C. 
This  application  of  this  sensitivity  is  currently  under  development. 

5  Summary 

A  new  adjoint  code  has  been  developed  for  the  NHL  Atmospheric  Variational  Data 
Assimilation  System  (NAVDAS),  which  performs  a  three-dimensional  variational  at¬ 
mospheric  analysis  in  observation  space.  The  NAVDAS  adjoint  can  be  used  to  effi¬ 
ciently  calculate  the  gradient  of  a  forecast  costfunction  with  respect  to  the  complete 
vector  of  observations  used  for  assimilation.  In  addition,  sensitivity  to  the  background 
and  the  specified  observation  error  can  be  calculated.  For  the  dry  observed  variables 
(temperature,  wind,  height),  the  observation  sensitivity  is  shown  to  be  accurate  to 
within  about  5%  of  the  theoretical  best-possible  result.  The  error  in  this  calculation 
is  likely  due  to  an  inconsistency  between  the  sensitivity  to  the  NOGAPS  initial  sur¬ 
face  pressure  field  and  the  NAVDAS  analysis,  which  provides  temperature  and  wind, 
but  not  an  analyzed  surface  or  sea-level  pressure.  At  present,  the  NAVDAS  adjoint 
does  not  include  sensitivity  to  moisture  observations,  but  it  is  configured  for  radiance 
assimilation  (e.g.,  sensitivity  to  brightness  temperature). 

This  report  describes  the  main  features  of  the  NAVDAS  adjoint,  including  key  sub¬ 
routines  and  the  sequence  of  steps  required  to  calculate  observation  sensitivity.  The 
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main  difference  between  the  forward  and  adjoint  versions  of  NAVDAS  is  the  trans¬ 
position  of  the  post-multiplier  P5H'^  from  the  forward  NAVDAS  into  the  operator 
HPft,  which  is  used  at  the  beginning  of  the  NAVDAS  adjoint  sensitivity  calculations. 

The  input  for  the  NAVDAS  adjoint  is  an  analysis  sensitivity  vector,  provided  by 
the  adjoint  of  a  forecast  model  such  as  NOGAPS.  A  special  interpolation  procedure  is 
required  to  transfer  the  analysis  sensitivity  from  the  NOGAPS  grid  onto  the  NAVDAS 
analysis  grid.  The  NAVDAS  adjoint  uses  the  imioA^tion  file  and  the  rsd_vec  file 
produced  by  the  regular  (forward)  NAVDAS  procedure. 

In  Section  4,  four  examples  of  observation  sensitivity  applications  axe  described: 
(i)  simple  examination  of  observation  sensitivity  gradients,  (ii)  impact  of  innovations 
on  actual  forecast  error  differences,  (iii)  sensitivity  to  h5q)othetical  observations,  and 
(iv)  targeted  observing.  These  examples  illustrate  the  capability  of  the  NAVDAS 
adjoint  to  provide  sensitivity  information  for  predictabihty  studies  and  to  study  the 
performance  of  the  data  assimilation  and  quality  control  procedures. 
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Appendix  A:  Scripts  and  Code 

Selected  scripts,  source  codes  and  makefiles  used  to  compile  and  run  the  NAVDAS 
adjoint  are  provided  on  the  NRL  computer  hadley  in  the  directory: 

/users_hadley/langland/NAVDAS_adjoint/.  The  NAVDAS  adjoint  currently 
runs  on  the  SGI  AMS  cluster  maintained  by  Fleet  Numerical  Meteorology  and  Oceanog¬ 
raphy  Center.  There  are  additional  input  and  data  look-up  files  located  on  AMS 
that  are  required  to  run  the  NAVDAS  adjoint, 
makefile 

NAVDAS_adj.f  (source  code) 
run_navdas_adjoint  (c-shell  script) 

NAVDAS_adj.ksh  (k-shell  script) 
gradtp.f  (interpolation  routine) 

track.dl-f-tl_GrV_TRACK  (targeted  dropsonde  input  file) 

For  additional  information  contact  Dr.  RolfLanglandat:  langland@nrlmry.navy.inil. 
Appendix  B:  Derivation  of  Observation  Impact  Equation 
Consider  two  nonlinear  forecasts  of  lengths  /  and  g  both  verifying  at  a  time,  t, 
for  which  there  is  a  verifying  analysis,  Xj.  Forecast  xy  is  started  from  an  analysis 
Xq,  and  forecast  x^  is  started  6-h  earlier  (corresponding  to  the  6-h  interval  of  the 
NAVDAS  data  assimilation  cycle).  The  6-h  forecast  of  x^  provides  the  first-guess  (or 
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background,  xj)  for  the  analysis  Xa  that  represents  the  initial  conditions  for  x/.  The 
difference  between  the  errors  of  forecasts  Xf  and  Xg  is  given  by  Ae^  =  Sf  —  eg.  We 
wish  to  derive  expressions  (or  ’’impact  fimctions”)  that  can  be  used  to  estimate  Ae^ 
using  adjoint  sensitivity  gradients  in  both  the  gridpoint  space  of  the  forecast  model 
and  the  observation  space  of  the  data  assimilation  procedure. 

We  first  define  quadratic  measures  of  the  two  forecast  errors  using  the  expressions 

e/  =  (  (x/  -  xt)  ,  C  (x/  -  xt))  ,  (Bl) 

^9  =  <  (Xfl  -  X() ,  C  (Xj  -  Xf))  ,  (B2) 

where  C  is  a  matrix  of  energy  weighting  coefficients  for  dry  total  energy  (Rosmond 
1997).  Using  (Bl)  and  (B2)  we  define  two  costfunctions 

=  Xf)  ,  C  (x/  -  Xf))  ,  (B3) 

—  (xj  —  Xf) ,  C  {xg  —  Xf))  ,  (B4) 

and  the  corresponding  first  derivatives 


||  =  C(x,-x,), 

(B5) 

g=C{x,-x.). 

(B6) 

The  true  value  of  e/  -  eg,  which  we  will  call  Ae^,  is  defined  as  (Bl)  -  (B2).  We 
wish  to  have  an  expression  for  Ae®  that  involves  sensitivity  gradients,  and  this  can 
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be  written  using  the  gradients  defined  in  (B5)  and  (B6)  as 

+  .  (B7) 

It  is  easily  shown  that  the  RHS  of  (B7)  is  equivalent  to  the  difi'erence  in  forecast 
errors  ef  —  Cg.  Note  that  it  is  necessary  to  use  sensitivity  gradients  involving  both 
the  /  and  g  trajectories  because  c/  —  is  the  difference  between  two  quadratic  error 
expressions. 

The  difference  between  the  forecast  trajectories  /  and  g  at  initial  time  is  Xq  —  x^, 
the  analysis  minus  the  background.  If  we  assume  that  this  initial  difference  evolves  in  a 
tangent  linear  sense  into  a  good  approximation  of  the  forecast  difference  Xf  —  Xg,  then 
we  can  estimate  e/  —  eg  using  Xq  —  and  sensitivity  gradients  which  the  forecast 
model  adjoint  has  mapped  back  to  initial  time  along  the  two  forecast  trajectories. 
Thus,  the  adjoint  model  maps  dJf  /  dxf  into  dJf  /  dxa  along  trajectory  /  and 
dJg  /  dxg  is  mapped  into  dJg  /  dxb  along  trajectory  g,  which  allows  us  to  write 

^e?=((x,-x.).g  +  g).  m 

Equation  (B8)  provides  an  estimate  of  e/  —  in  the  gridpoint  space  of  the  fore¬ 
cast  model.  The  result  is  not  exact,  although  Ae®  in  (B7)  is  exact,  because  the 
sensitivity  gradients  in  (B8)  are  approximations  obtained  with  the  adjoint  of  NO¬ 
GAPS.  Prom  (B8)  we  can  derive  an  equivalent  expression  for  in  the  observation 
space  of  the  analysis  procedure.  First,  we  recognize  from  Eq.  (6)  in  Section  2.2 
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that  Xa  “  X5  is  equivalent  to  K  (y  —  Hx^,),  where  the  Kalman  gain  matrix  K  = 

+  R]  (Daley  and  Barker,  2001).  We  may  therefore  re-write  (B8) 
as 

Je?  =(K(y-HxJ,g  +  g).  (B9) 

Now,  since  K  is  linear,  we  can  use  the  general  definition  of  an  adjoint  operator 
(Ky  ,  x)  =  (y  ,  K^x)  to  write 

+  (BIO) 

The  operator  represents  the  adjoint  of  the  data  assimilation  procedure,  which 
defines  a  sensitivity  gradient  in  observation  space 


And  using  (Bll)  to  substitute  into  (BIO)  we  obtain  (16)  used  in  Section  4.2 

/  dJ^\ 

M  =  ((y-Hx6)  ,  ,  (B12) 

which  involves  only  observation  space  quantities.  If  there  are  no  observations, 
then  Xq  =  X(,  and  (5e®  =  0.  The  accuracy  of  the  observation  impact  calculated  in 
observation  space  (B12)  is  essentially  equivalent  to  the  gridpoint  space  calculation 
(B8). 
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Appendix  C:  Derivation  of  Sensitivity  to  Observation  Error 

The  operator  representing  the  solver  step  in  the  forward  NAVDAS  is 

M  =  [HPf.H'^  +  R]  .  (Cl) 

Using  M  the  ”cob”  vector  (  z  )  is  defined  as 

z  =  M-^  d  ,  (C2) 

where  d  represents  the  innovation  vector  (y  —  Hxb) .  Next,  (C2)  is  re-witten  as 

M  z  =  d  ,  (C3) 

and  in  perturbation  form  as 

<JM  z  +  M  (iz  =  5d  .  (C4) 

In  (C4)  6z  represents  perturbations  of  the  ”cob”  vector,  and  5M  represents  per- 
tmbations  of  the  M  operator,  which  could  include  changes  to  the  specified  observation 
or  background  error  covariance.  Next,  (C4)  can  be  re-written  as 

M  (Jz  =  (5d  -  (JM  z  ,  (C5) 


and. 


Sz  =  M-^  ((Jd  -  5M  z) 


(C6) 


FVom  (C6),  and  noting  that  M  is  self-adjoint,  we  obtain  the  following  adjoint 
equations 
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where  (*)  denotes  a  vector  Shur  product.  In  (C7)  the  innovation  sensitivity  dJ  /  dd 
is  equivalent  to  {dJ  /  dy)  if  the  background  is  not  perturbed.  We  can  thus  substitute 
from  (C7)  into  (C8)  and  obtain 


dJ  _  ^ 
dM~  dy 


(C9) 


The  quantity  dJ  /  5M  is  a  vector  in  observation  space  that  represents  the  sen¬ 
sitivity  of  J  with  respect  to  the  operator  M.  Here,  however,  we  are  just  concerned 
with  sensitivity  to  the  elements  Tj  of  the  observation  covariance  matrix  R.  Then, 
assuming  the  operator  Pj,  is  not  perturbed  and  since  R  is  diagonal  we  obtain 


dJ  dJ 

dr,- 


(CIO) 


which  is  (23)  in  Section  4.5. 


